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Abstract. We present a simple and easy to implement method for the numer- 
ical solution of a rather general class of Hamilton- Jacobi-Bellman (HJB) equa- 
tions. In many cases, classical finite difference discretisations can be shown to 
converge to the unique viscosity solutions of the considered problems. How- 
ever, especially when using fully implicit time stepping schemes with their 
desirable stability properties, one is still faced with the considerable task of 
solving the resulting nonlinear discrete system. In this paper, we introduce a 
penalty method which approximates the nonlinear discrete system to an or- 
der of 0(l/p), where p > is the penalty parameter, and we show that an 
iterative scheme can be used to solve the penalised discrete problem in finitely 
many steps. We include a number of examples from mathematical finance for 
which the described approach yields a rigorous numerical scheme and present 
numerical results. 
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Introduction 

In their most general form, Hamilton-Jacobi-Bellman (HJB) equations arise when 
applying Bellman's principle of optimality to problems of stochastic optimal control, 
in which case the solution to the HJB equation can often be found to solve the 
control problem (cf. }26|). There is a wide range of practical applications to optimal 
control, amongst others in mathematical finance, and many examples are listed in 
[91 HH] . Naturally, the high practical relevance gives rise to the question of how 
to best solve stochastic control problems numerically. Early (and still relevant) 
methods rely on Markov chain approximation, of which an extensive overview can 
be found in jTTl [TBI [9] . The notion of viscosity solution, introduced in [14] , together 
with the most helpful results of [2] , permits another line of attack, namely the direct 
solution of the corresponding HJB equations by discretisation methods; for details, 
see [2J HI |3] and references therein. It is interesting to note that both approaches 
can be found to be closely related, which is most apparent for explicit time stepping 
methods (e.g. cf. [23]). 

In this paper, we are concerned with the numerical solution of a certain kind of 
HJB equation with a finite control set. Using results of [2J, where it was - roughly 
speaking - shown that every stable and consistent discretisation converges to the 
financially relevant solution, one can often find numerical approximations by em- 
ploying standard finite difference schemes. However, the resulting discrete system 
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Framework considered 
in Section 1 



Continuous HJB Equation 





Convergence to the unique viscosity 
solution. 



Modified 
Newton 
Scheme 




Discretised HJB Equation 



Approximate to an order of 0(l/p), 
where p>0 large. 



Penalty Approximation 
of the Discretised HJB 
Equation 



Figure 1. To numerically approximate the HJB equation in (A), we 
consider a discrete version (B). Given some conditions that have to be 
checked for every particular example, the grid convergence of (B) to (A) 
can be guaranteed. We then approximate (B) by the penalty problem 
(C), which in turn can be solved exactly using algorithm (D). 



is only straightforward to solve if fully explicit time stepping is used, in which 
case one suffers from undesirable stability constraints. When using fully implicit 
or weighted time stepping schemes, one faces a nonlinear discrete system, which 
is generally nontrivial to solve. Existing algorithms are, in one way or another, 
all related to policy iteration techniques, which have been used in the setting of 
Markov chain approximation [THJ, [H] as well as for HJB equations [TOl 113 US] • In 
the present paper, we introduce a novel approach by which the nonlinear discrete 
system can be solved using a penalisation technique. More precisely, we show how 
a penalty approach that was devised for American option problems in |llj can be 
extended to solve discretised HJB equations provided the discretisation matrices 
obey a certain structure. (Penalty methods in general have a long history and 
have, for example, already been used extensively in the setting of variational in- 
equalities in [5].) Examples from mathematical finance that fit into our framework 
include uncertain volatility models [20], transaction cost models [19], and unequal 
borrowing/lending rates and stock borrowing fees [HE]; a succinct overview of these 
and similar models can be found in [101 112| . 

Structure of this Paper. We aim to devise a penalty-based numerical scheme 
for the solution of a rather general class of HJB equations; the broad concept is 
depicted in Figure [T] 

Section [7J We define the considered class of HJB equations (denoted by (A) in 
Fig.[lJ and - having fully implicit and weighted time stepping schemes in mind - a 
fairly general nonlinear discrete problem (denoted by (B) in Fig.[T]). The structure 
of (B) is chosen such that we have an unconditionally stable numerical scheme 
which converges to (A) if the latter has a unique viscosity solution and satisfies a 
strong comparison principle. 

Section^ The nonlinear discrete problem (B) is nontrivial to solve and will be our 
main consideration. In a first step, we introduce a penalised problem (denoted by 
(C) in Fig. [I]) which approximates (B) to an order of 0(1/ p), where p > is the 
penalty parameter. It is worth pointing out that as the discretisation in (B) is stable 
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(due to the monotonicity of the discretisation matrices), any errors introduced by 
the penalisation will stay within bounds and, hence, will not compromise the overall 
numerical scheme. Since (C) is a nonlinear discrete problem itself, the feasibility 
of the penalty approach crucially depends on the possibility of conveniently solving 
(C) in a second step; in Section [3j we present an algorithm doing exactly that. 

SectaonfJ} We present an iterative procedure (denoted by (D) in Fig. [I]) with finite 
termination which finds the exact solution to ( C). Even though the theoretical 
bound on the number of iterations depends on the the number of grid points and 
the cardinality of the finite control, the algorithm always converges - independently 
of the grid size and the penalty parameter - in less than four steps for our numerical 
examples. 

Section [7} We give a detailed outline of how various examples from mathematical 
finance fit into our framework. Furthermore, for the HJB equation arising in Eu- 
ropean option pricing with different borrowing/lending rates and stock borrowing 
fees, we show in detail how convergence of (B) to the unique viscosity solution of 
(A) can be guaranteed. 

Section [5| We conclude by presenting numerical results and compare our approach 
to the method of policy iteration. Also, we provide a brief proof showing that 
policy iteration, just like the iterative scheme proposed by us in Section |3j has 
finite termination when applied to the kind of HJB equation considered in this 
paper. 

1. Problem Formulation 

We begin by (slightly informally) introducing the type of HJB equation (see Prob- 
lcm |l.l[ ) considered in this paper and for which we subsequently develop a fully 
coherent numerical scheme. In Section [4j we will then show that the problems from 
mathematical finance listed in the introduction all match our framework and that 
rigorous problem formulations can be obtained by employing the theory of viscosity 
solutions (as introduced in [H]). 

Problem 1.1. Let §, |S| < oo, be a finite set of controls. Let f2 C ffi" (with n e N) 
be an open set. Let C s , s G §, be a number of linear differential operators on fi. 
Find a function V : Q — > R such that 

(1.1) mm{C s V : s G §} = 0. 

At this stage, we are not yet concerned with possible boundary/initial conditions 
needed for the uniqueness of a solution to Problem |1.1| since formulation (1.1 1 is 



primarily motivational. Nevertheless, we point out that (if we have an additional 
linear differential operator £*) formulation ( |1.1| is sufficiently general to include 
problems of the form 

(1.2) C*V + xmn{C s V : s e §} = 



when defining C s := C s + C* , s£§. In many applications, C*V in (1.2 1 is taken 



to be the time derivative of V (i.e. C*V = dV/dt). Also, the approach described 



in this paper works analogously when replacing the "min" in equations (1.1 1 and 



(1.2) by "max", which will be explained in Remark 1.4 



Remark 1.2. Let N E N and define Af := {1, . . . ,N}. Consider two vectors x, 
y G R N and a matrix A G R NxN . For any i £ Af, we denote by (x)i and (A)i the i- 
th coordinate and the i-th row of vector x and matrix A, respectively. When writing 
x > 0, we generally mean that (x)i > for all i G Af, and we take z :— min{x, y} to 
be the vector satisfying {z)i — min{(ai)j , (y)i} for all i G Af. The definitions extend 
trivially to other relational operators and to the maximum of two vectors. 
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In the following, whenever the choice of N £ N is clear from the context, we will 
use the set AT as introduced in Remark |1.2| without explicitly saying so. 



We make the assumption that we can find a sensible (i.e. stable and convergent) 
discretisation of Problem |1.1| that results in discrete systems of the following form. 

Problem 1.3. Find x £ R N such that 

(1.3) mm{A s x - b s : s € §} = 0, 

where 

• b s £ R , s£§, are vectors, 

• A s £ R NxN , s £ S, are matrices with non-positive off-diagonal entries, 
positive diagonal entries, non-negative row sums, and at least one positive 
row sum, 

• and the positive row sum is equally located in all matrices A s , s£§. 

In Section |4.1| we will show that stable and convergent discretisations satisfying 



the setup of Problem 1.3 can be found in many relevant cases (and, most notably, 
for many examples from finance), justifying our assumptions. The conditions on 
the matrices are motivated by the fact that a square matrix with non-positive off- 
diagonal entries, positive diagonal entries, non-negative row sums, and at least one 
positive row sum is an M-matrix (cf. [U I24j). which means, in particular, that it 
has a non- negative inverse; we shall make frequent use of this fact. 



Remark 1.4. If the "min" in equation (1.1) is replaced by "max", equation (1.3) 
becomes 

(1.4) imn{A s x - b s : s e §} = 0, 

with A s := —A s and b s := —b s . In this case, —A s , s £ S, is an M-matrix, which 
means that A s has non-positive inverse A -1 < 0; therefore, all proofs that follow 
(Theorem 1.5 and Sections^ and[3|) can be shown to hold for both situations as 



soon as they hold for one. However, for the sake of simplicity, we solely focus on 



situation ( 1.3 ) 



Theorem 1.5. There exists a unique solution to Problem \l.S[ 



Proof. The proof of existence will be given in Corollary |2.5| Suppose we have two 
solutions x\ and x 2 ■ For every i £ N ', there exists an Sj £ S such that 

{A Sl ) l x 2 - (b Si )i = 0, 

and it also is 

(A s A lXl - (b Si )i > 0. 
Denote by A* £ R NxN the matrix consisting of the rows (A Si )i , i £ AT, which is 



an M-matrix by our assumptions on the discretisation (cf. Problem 1.3 ). We have 

A*( Xl -x 2 ) > 0, 

and x\ — X2 > follows since (A*)^ 1 > 0; swapping the arguments, we get x 2 — x\ > 
0, which then proves the theorem. □ 

2. Penalisation of the Discrete Problem 

In this section, we discuss how the nonlinear discrete system in Problem |1.3| can be 
approximated by another (also nonlinear) problem using a penalisation technique, 
and it will turn out that the penalised formulation can very efficiently be solved 
numerically. Our penalty approximation to the discrete HJB equation can be seen 
as an extension of the penalty scheme for American options as introduced in [llj . 
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2.1. The Penalised Problem. Our penalty approximation looks as follows. 
Problem 2.1. Let so G S, p > 0. find .x p G R N such that 

(2.1) (4,„ &„-&,„)-/> J2 max{& s - A, x p , 0} = 0. 

ses\{s } 



In this, p > is called the penalty parameter, and ( 2.1 ) can be thought of as apply- 
ing a penalty of intensity p whenever an expression A s — b s x p , s £ §\{so}; turns 
negative; by pushing up the expression p X) s gs\{ s „} max{6 s — A s x p , 0}, the penalty 
parameter causes (A SQ x p — b S() ) and, implicitly, x p to grow until no constraints are 
violated anymore. 

Theorem 2.2. There exists a unique solution to Problem\2.1\ 



Proof. The proof of existence will be given in Theorem |3.5| Suppose we have two 
solutions x P: i and x Pt2 ■ It is 

(2.2) (A sg x Pt i - b sg ) - p max{& s - A s x pA , 0} = 

ses\{s } 

and 

(2.3) (A So x p , 2 -b So )- p max iA - A s x p ,2 , 0} = 0, 

ses\{ So } 



and, subtracting (2.3 1 from (12. 21), we get 



(2.4) A So (x Pi i - x p , 2 ) - p Y max{& s - A s x pA ,0} 

ses\{s } 

+ p max{& s - A s x p ,2 , 0} = 0. 

s£S\{s } 

Now, let i G M and r G S\{so}, and consider 

(2.5) ( — pmax{fe r — A r x Pt \ , 0} + pmax{& r — A r x p ^ 2 , 0}) . 

= -pmax{(6 r )i- (A r x pA )i ,0} + pma,x{(b r ) l - (A r x p , 2 )i,0}, 
where we distinguish between the following different cases. 

(i) If max{(b r )i — (A r x p ,i)i , 0} = max{(6 r )i — (A r x Pj2 )i , 0} = 0, then expres- 



sion (2.5) is identical to zero, 
(ii) If max{(6 r )j — {A r x p,\)% j0} > and max{(6 r ) i — (A r x Pj2 )i , 0} > 0, then 
expression ( |2.5| can be written as p(A r x p ,i)i — p(A r x Pl2 )i = (pA r (x Pl i — 

(hi) If max{(6 r )i — {A r x Py i)i ,0} > and max{(6 r )i — (A T x Pl2 )i , 0} = 0, then 



expression (2.5) is negative, i.e. —p{b r )i + p(A r x Pi i)i < 0. 



(iv) If max{(6 r )j — (A r x Pt i)i >0} = and max{(b r )j — (A r x Pt2 )i , 0} > 0, then 



expression (2.5) equals 

p(b r )i - p(Ar X Pt2 )i 

= p(W)i - p(A r x Pt2 )i + p(A r x Pt i)i - p{b r )i - (p(A r x p ,i)i - p{b r )i) 

= p{A r {x pA - x Pi2 )). + p(b r )i - p(A r x Ptl )i , 

in which p(b r )i — p(A r x Pt i)i < 0. 
Now, noting that the M-matrix property of A So is preserved if we add rows of the 



matrices A s , s £ S\{so}> to some of its rows, expression ( 2.4 1 implies, in the light of 
the four cases discussed above, that there exists an M-matrix A* G ~R NxN satisfying 

A*(x Pt i - x p ,2) > 0, 
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which then gives (x Pt i — x Pt 2) > since (A*) 1 > 0; by the same reasoning, we can 
also get (x Py 2 — %p,i) > 0, which then completes the proof. □ 

2.2. Properties of the Penalised Problem. What follows are properties of the 
penalised problem. We prove that the penalty approximation converges to the 
original problem as p — >• oo and provide an estimate for the penalisation error. 



Lemma 2.3. Let x p be the solution to Problem 2.1 There exists a constant C > 0, 
independent of p, such that ||xp||oo < C. 

Proof. We rewrite equation ( |2.1[ ) to get 

(2.6) A So x p = p max{b s - A s x p , 0} + b So . 

s€S\{s Q } 

Hence, for every component i £ Af, we have that 

(A So x p )i = (b So )i or 3 s £ S\{s } s.t. (A s x p )i < {b s )i . 

Therefore, we may, for i £ Af, define £ S to be such that (A Si x p )i < (b Si )i is 
satisfied, and, introducing A* £ M. NxN to be the matrix having as i-th row the i-th 
row of A s . , i £ Af, and defining b* £ M. correspondingly, it is 

A* x p < b*. 

Since A* is an M-matrix, it is (A*) -1 > 0, and we get x p < (A*)^ 1 ^, in which 
(A*)^ 1 and b* can be bounded independently of p because there are only finitely 
many compositions that can be assumed. To see that the negative part of x p is 
bounded as well, we note that from (12. 61) we have 



Xp > — (A So ) 1 b S(l , 

in which (A^) -1 > 0. □ 



Lemma 2.4. Let x p be the solution to Problem 2.1 There exists a constant C > 0, 
independent of p, such that 

(2.7) A S0 x p - b S0 > 0, 

(2.8) A s xp - b s > - -, s£ S\{s }, 
and, for every i £ Af, it is 

(2.9) (A^Xp-bsJ^O 
or 3 s £ S\{sq} such that 

(2.10) {AsXp-bs), < -. 



Combined, properties ( 2.7 1 — ( 2. 10 1 imply that \\ min{^4 s x p — b s : s £ §}||oo < C / p. 
Proof. For s £ §\{so}, it is 

p(b s — A s x p ) < p max{6, r — A r x p , 0} 

r£S\{s } 

— A So Xp b S() , 
and, applying Lemma |2.3| to the last expression, we get that 

(2.11) A s x p -b s >-- 

P 

for some constant C > independent of p. Furthermore, we have that, for every 
component i G AT, it is 

(A Sg x p — b So )i = 
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or 3 s e §\{s } such that 

(b s - A s x p )i > 0, 
which, together with (2.111, extends to 

|(^x-6.)<| < -. 



□ 



Corollary 2.5. Let x p be the solution to Problem 2.1 As p — » oo, x p converges to 
a limit x* which solves Problem \l.S\ 



Proof. Since (£ p ) p >o is bounded (as seen in Lemma 2.3 1, it has a convergent sub- 
sequence, which we will not distinguish notationally; we denote the limit of this 
subsequence by x* . From (2.7) and (2.8), we may deduce that 

A s x* -b s > 0, s e S. 



Now, consider a component i € AT. From (2.9) and (2.10), we know that, for every 
p > 0, there exists an s^ p € S such that 



\(A Si p x p - b Si Ji\ < 



C 
9 ' 



Recalling that |S| < oo, we can then infer that there exists once more a subsequence 
of (x p )p>o , again not notationally distinguished, and an s* £ § such that 

(2.12) \(A s ,x p -b st )i\<- V P >0, 
which means 

(2.13) \(A S * x*-b st )i\ =0. 

Altogether, we see that x* solves Problem |1.3| Since Problem |1.3| has a unique 



solution (cf. Theorem 1.5 ), the just given result holds not only for subsequences of 
{Xp)p>o * but the whole sequence converges. □ 

Theorem 2.6. Letx p and x* be the solutions to Problems 2.1 and \1.3\ respectively. 
There exists a constant C > 0, independent of p, such that 

C 



< 



P 



Proof. We continue in the setup of the proof of Corollary 2.5 From (2.12) and 
(2.13), we know that the original sequence (a; p ) p >o of solutions to Problem 2.1 has 
a subsequence (s p ') P '>o such that, for every component i £ J\f, there exists an 
s* e § satisfying 

\(A„*)i (x p > -x*)i\ < ^ Vp' >0. 

Denote the matrix consisting of the rows {A s * )j , i E Af, by A* 6 M. NxN . Noting 
that A* is an M-matrix and hence invertible, we obtain that 



(2.14) 



< 



K^niooC 



v P ' > o, 



in which || (A*) -1 can be bounded independently of p because there are only 
finitely many compositions it can assume. Since, by the same arguments, every 
subsequence of the original sequence {x p ) p> o can be shown to have a subsequence 
satisfying (2.14), we deduce that the whole sequence satisfies 
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for some constant C > independent of p. □ 



3. An Iterative Method for Solving the Penalised Problem 

Having seen in the previous section how penalisation can be used to approximate 
Problem |1.3[ it remains to answer the question how to numerically solve the non- 



linear discrete system (2.1 1 appearing in the penalised formulation. To this end, we 
present an iterative method - distantly resembling Newton's method - which has 
finite termination and converges to the solution to Problem |2.1| In the situation of 
American options, a similar algorithm was used in |llj . 

We are trying to find y € which satisfies ( |2.1| ; hence, defining 

G(y):={A So y-b So )-p £ max{6 s - A s y, 0}, y e R N , 

s<£S\{s } 

we need to solve G(y) = 0. 

Remark 3.1. Let £ be a vector (or a matrix) in R N (or M. NxN ). For y £ M. N , we 
denote by £f , s£ §\{so}, the vector (or matrix) consisting of 

• all rows of £ where the corresponding rows of b s — A s y are positive 

• and where all other rows are set to zero. 

Algorithm 3.2. Define 

j G (y):=A S0 + P AV ^ y^ RN - 

s£S\{s } 

Let x° £ R N be some starting value. Then, for known x n , n > 0, find x n+1 such 
that 

(3.1) J G (x n )(x n+1 - x n ) = -G(x n ). 

Clearly, the just introduced algorithm bears strong resemblance to Newton's method. 
However, since (as we will see) it converges globally in finitely many steps, it is 
more powerful than standard Newton iteration, for which convergence can only be 
guaranteed in a sufficiently small neighbourhood of the solution. Furthermore, the 
function G(-) is not differcntiablc, which means standard Newton iteration is not 
directly applicable. For general results on extensions of Newton methods to semi- 
smooth Newton methods for the solution of particular non-smooth equations, see 
[131 [TB] and references therein. 

The algorithm is well defined since Jq is a non-singular matrix. Also, equation 
d3. Ih can be written as 



(A S0 +p Af)(x n+1 ~x n ) = -(A So x n -b So ) + p (€~Afx n ), 

s&\{s } s£S\{s } 

which is equivalent to 

(3.2) (A S0+P J2 Af)x n+1 =b SB+P J2 b f- 

s<£S\{s } sSS\{so} 

Lemma 3.3. Let x° € R N be some starting value, and let (x Tl )^ =0 be the sequence 



generated by Algorithm 3.2. Lt is x n < x n+1 for n > 1 
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Proof. Writing (|3.2|) for x n and x n+1 , we obtain 

b. 



n+l 



(a so + p y A * 

ses\{s } 

and {A S0 +p Y Af' 1 )x n = b s 
ses\{s„} 

and subtracting yields 

(3.3) (A S0+P J2 Af)(x n+1 -x n ) 
ses\{s„} 

= P y b f-p E A T xn -p E bX s 

s6S\{s } sGS\{s } sGS\{s } 



p E b * 

s£§\{s } 

p E b * 

sGS\{s } 



P y A f~^ x " 

sG S\{s } 



In this, A So + pX)ses\{ So } i s an M-matrix (and has a non-negative inverse), 
and, therefore, the proof is complete if we can show that the right hand side of 
expression (3.3) is non-negative; to do this, we consider the rows i G Af of the right 
hand side of ( |3.3| separately. Let s £ §\{s }. 

(i) If {Af )i = and (Af = 0, it is p{bf )< - p(Af x"), - p(&f + 
p(Af V); = 0. 

(ii) If (Af )i > and (Af > 0, it is (Af ), = (Af and (hf )< = 
(&f , and, hence, it is p{bf ) l -p{Af x n ), -p{bf ' 1 ) l +p{Af ~\ n )i = 0. 

(hi) If {Af)i > and (Af = 0, it is p{bf) t - p(Af x n \ - p{bf~\ + 

p{Af- 1 x n ) i = p(6f )i - p{Afx% > 0. 
(iv) If {Af)i = and {Af^i > 0, it is p{bf)i - p{Af x n \ - p{bf \ + 

p{Af x n )i = —p{bf )i + p(Af x n ) t > 0, where the last inequality 

holds because max{(6 s )i — {A s )i x n , 0} = 0. 

This completes the proof. □ 

Lemma 3.4. There exists a number neN such that, for every starting value x° , 
Algorithm\3.2\ terminates after at most k steps. 



Proof. Since A So + P^2 se s\{ So } A x s is a non-singular matrix, we may write (3.2 ) as 
(3.4) x n+1 = {A S0+P Y, Af)-\b S0 +p Y b f)- 

s£S\{s } sGS\{s } 



There is only a finite number of compositions for the right hand side of (3.4) as it 
consists of given matrices of which some rows are set to zero. Hence, there exists a 
set C, \C\ < oo, C independent of x° , such that x n € C for all n G N. Also, by the 
properties of the algorithm, if x n +1 = x n for some n' £ N, then it is x n — x n for 
all n > n' . This, together with the monotonicity of the algorithm (cf. Lemma 3.3 1, 
proves the lemma. □ 



Theorem 3.5. Independently of the starting value x° , Algorithm 3.2 converges (in 



finitely many steps) to a limit x* which solves Problem 2.1 



Proof. As the proof of convergence has already been given in Lemma 3.4 it only 
remains to show that the limit x* satisfies G(x*) = 0. This follows easily from 
the fact that x* — x n+1 = x n for all n > k, which means that (3.1) reduces to 
-G{x*) = ~G{x n ) = J G {x n )(x n+1 - x 1 



0. 



□ 
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4. Examples of HJB Equations Arising in Finance 

Let r, q > 0, a > denote interest rate, (continuously paid) dividend rate and 
volatility, respectively. We introduce the Black-Scholes operator 

(4-1) C& V := ^ + -*S — 2 + (r - q)S- - rV, 

where V = V(S, t) is some function of S and t. 



Remark 4.1. Let M, N G N. Suppose we want to solve equation (4.1) numerically 
backwards in time on the interval [0,S max ] x [0,2^]. Also, suppose that a terminal 
value and boundary conditions have been fixed. Let time and space grids be given 
by {jk : j + 1 G A4} and {ih : i + 1 € M}, respectively, where A4 := {1, . . . , M}, 
Af := {1, . . . , N}, k := T / [M — 1) and h := S max /(N — 1). Using a fully implicit 
finite difference method with a one-sided difference in time and central differences 
in space, we obtain the following scheme: For every time step j G A4\{1}, we have 
to find V^ 1 G R N such that 

,i; ; ri v ; + V j = 0, 

where G M. N is the vector known from the previous time step and A T g'g' 7 G R JVxAr 
is a tridiagonal matrix with (bi)(L 1 on the diagonal, (ci)^^ 1 on the upper diagonal, 
and {ai)^ =2 on the lower diagonal; more precisely, for i G A/"\{1, N}, the coefficients 
are given by 

a r t ' q ' a = -^i 2 a 2 k+^i(r - q)k, 



h r,q,a 



= 1 + i 2 <j 2 k + rk, 



and c i ' q ' (J = —-i 2 o~ 2 k — -i(r — q)k, 

and b r { q,a , c[' q ' <T , a r ^ ,a and are chosen to accommodate the boundary condi- 

tions. In finance, boundary conditions can usually be found such that A^g 7 is an 
M-matrix. 

Proof. See [35] ■ □ 

The examples from mathematical finance listed in Section [l] (uncertain volatility 
models [20], transaction cost models [19], and unequal borrowing/lending rates and 



stock borrowing fees [TJ [HI H31)j ah result in HJB equations of the form (1.1) where 



every differential operator C s , s G S, can be interpreted to be of Black-Scholes type 



(4.1 ) for some set of parameters r, q and a. (For details, we refer to [UIIIH], where 



the models have been described at length.) As seen in Remark 4.1 a finite difference 



discretisation of (4.1) results, for every timestep, in a linear system of equations 
where the defining matrix is an M-matrix, and, therefore, the above problems all 
fit into the framework of Problem |1.3[ a detailed outline of one representative case 
will be given in the next section. 

4.1. Example: European Option Pricing with Stock Borrowing Fees and 
Unequal Borrowing/Lending Rates. We would like to price a European op- 
tion that pays P(S) at time T > 0, where S denotes the price of the underlying. 
We assume the standard Black-Scholes model with the extension that the cash 
borrowing rate (denoted by r^) and the cash lending rate (denoted by r{) are not 
necessarily equal, and, also, that there are stock borrowing fees of r/ reflecting the 
fact that shorting the stock may have a certain cost. To preclude arbitrage, it must 
be r\j > ri > r/ . Different borrowing/lending rates have been discussed in [TJ[5], 
and stock borrowing fees are studied in [T3] ; the resulting mathematical model has 
been explained in |10j . 
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To price a short position in the option, we have to solve the following nonlinear 
problem. We define S max > to be a large upper limit on the value of the un- 
derlying S and set f2 := (0, S max ). For simplicity of presentation, we assume that 
P(0) = P(S max ) = 0; most other boundary conditions can be treated similarly. 

Problem 4.2. Let a > and r b >n>r f >Q. Find V : H x [0, T] -> R such that 
V(S, t) = P(S) for (S, (jefflx (0, T], V(S, T) = P(S) for SeH, and 

(4.2) max {C^/V : (r, q) e §>} = 

onQ x (0,T); here, it is § := {(r ; , 0), (r b , 0), (rj ,r f ), (r b , r b - n + r/)}. 

We use the notion of viscosity (sub-/super-) solution as developed in [T3], and, in 
that framework, the following strong comparison principle - which guarantees the 
uniqueness of any viscosity solution - can be shown to hold for the pricing problem. 

Theorem 4.3. Suppose u and v are, respectively, viscosity sub- and supersolution 
of Problem 4-2 If the payoff satisfies P G C(fl), then it is u < v on x (0, T]. 

Proof. After transforming the terminal value problem into an initial value problem, 
this can be found directly by checking the necessary conditions in [TJ] or [7]. □ 

It now seems natural to apply the discretisation technique from Remark |4.1| to every 
Black-Scholes operator in equation (4.2 1 individually to obtain a nonlinear discrete 
approximation of Problem |4.2| 



Problem 4.4. Let M , N £ N and use the notation of Remark 4-1 For a and § 



as in Problem \4-^ find {V^j^M C 1^ such that 
(4.3) mml^/V^'- 1 - V 3 : (r,q) G §} = 0, j G M\{1}, 

where 

• V M G K. is the function P evaluated on the grid {ih : i + 1 G TV} 

• and b r { q,G — b r I f' a = 1 and c^ q ' a = a^ 1 '' 7 = in all matrices A^g 7 . 

What follows are a few properties - stability, consistency and monotonicity, used 
as defined in [2] - of the just introduced numerical scheme; based on these, we will 
then be able to deduce convergence. 



Lemma 4.5. If a 2 > r b , the discretisation in Problem 4-4 * s unconditionally stable 
in the maximum norm. 

Proof. For any (r,q) G § and i G M\{1,N}, it is b r - q ' a > 1 and a[' 9 ' CT , c[' 9,CT < 0. 
Hence, from (4.3), for any given j G A^\{1} and i G N} there is choice of 

(r, q) G § such that 

\yJ~ 1 \ - 1 \yo _ cT'^V^ 1 - c r ' q ' a V j ^ 1 \ 

\ v i \— ,r,q,a\ v i u i v i-l c i * i+1 I 



< 1 Crg " max(|^|,|V^ 1 |,|V^ 1 1 |), 



and we infer that 



max \V?- X \ <max(|V A /" 1 |,|^" 1 |, max \Vf\)< max P(ih), j G M\{1}. 

Ki<N 1 Ki<N 1 ' Ki<N 



□ 



We point out that the constraint on the size of a required in Lemma 4.5 can be 
avoided when using upwinding instead of central differences for the space deriva- 
tives. However, since, in practice, the implicit scheme runs stably regardless of the 
choice of <r, we do not investigate this point any further. 



Lemma 4.6. The discretisation in Problem 4-4 * s consistent. 
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Proof. Let 4> : Vt x [0, T] -> R, and, for i e N and j G M, write <f>\ := <j){ih,jk). 
We then define 

&(h,k,ih,jk,^ ,<j>) := 

1^4-1 + V^ti + <"*"#fi ~ if * G TV} A j G M\{M}, 

\#'-P(tJ0 if*G{l,JV}vj = M, 

and, for general (5, i) E f2 x (0, T) - which is not necessarily on the grid - we set 

&(h,k,S,t,4>{S,t),4>) := &(h,k,ihjk,4>i ,</>) 

using i = \S/h~\h and j = \t/k~\k. This definition of the function & is in line with 
the notation in [2], and, clearly, our numerical scheme from Remark |4 . 1 1 corresponds 
to &(h, k,ih,jk, Vf , {V^j^M) = 0. Now, if <f> is smooth on Q, x [0,T], we have 

k 

0, t to and S — > So , where 

(£^V)(S ,i ) ii(S ,t )eflx(0,T), 

min{(£™^)(S ,t ) , t/>(S ,<o) - ^o)} if S G 90 or t = T, 

and the relation also holds when replacing "liminf", ">" and "min" by "limsup", 
"<" and "max" , respectively. This proves the lemma. □ 




Lemma 4.7. The discretisation in Problem 4-4 * s monotone 



Proof. We use the notation of Lemma |4.6| It is 

&(h,k, S,t, z,<f>) < &(h, k, S, t, z, tp) 
if <f> > tfi since a[' 9,<T , c i ,q,iy < 0. □ 

Finally, we have done enough preparatory work to formulate the following theorem 
on the convergence properties of our numerical scheme. 



Theorem 4.8. In the limit h, k —> 0, the solution (V^j^M to Problem 
verges uniformly on D, x [0, T] to the unique viscosity solution to Problem 



44 
4.2 



con- 



Proof. Using Lemmas |4.5| |4.6| and 4.7 this follows from results in [2]. □ 

5. Numerical Results 



In this section, we numerically solve Problem 4.2 using the penalty approach pre- 
sented in Section [2] To have a non-trivial financial instrument, we choose 

fS/4-25 if S G [100,200), 
P(S) = I -5/4 + 75 if Se [200,300), 
I else, 

which is the payoff function of a butterfly spread; unless stated differently, all other 
parameters are set as given in Table [T] 



n 


n 


r f 


a 


T 


Q 

^max 


M 


N 


tol 


P 


0.15 


0.1 


0.08 


0.4 


1 


600 


400 


400 


le-08 


le04 



Table 1. The parameters used in the numerical computations. (For 
illustrative purposes, we take unrealistically high interest rates.) 
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After discretising the problem (cf. Section 4.1 1, we are left with the nonlinear dis 



crete equation in Problem |4.4| (which corresponds directly to Problem 1.3 1. To 
solve it, we employ (i) the penalty approach from Section [2] (which means we solve 
Problem |2.1| instead) and (ii) the policy iteration as devised in [TU1 12] ; a detailed 
outline of policy iteration adapted to our case is given in Appendix[A] The following 



result is very similar to Lemma 3.4 



Proposition 5.1. There exists a number A € N such that, for every starting value 
x° , the method of policy iteration terminates after at most A steps. 

Proof. See Appendix [XJ □ 

Remark 5.2. The methods used to solve Problems and \2.1\ are of iterative 
nature, which means we need a valid check for convergence. 

• When solving Problem \2.1\ by Algorithm \3.2\ we abort the iteration as soon 
as we find an x™ such that 

\\(A S0 x n p -b S0 )- pJ2 se s\{s } max{6 B - A s x n p , O}^ < 
|| max{6 s : s G S}||oo 

• When solving Problem \l.S\ by policy iteration, we abort the iteration as soon 
as we find an x n such that 

(Ax n -b s )i > td e ^ x§ 
||max{6 s : s G S}|U " v ' 

and, for every i G N , there exists an r G § such that 
(A r x n b r )i 



max{6 s : s G §}| 



< tol. 



In Figures [2] and [3j we see the solution to Problem |4.2| for different values of r& , 
ri and rj. To understand the effects, it helps to see the butterfly spread as having 
a call side (the left hand side) and a put side (the right hand side). For a short 
position in a vanilla call option, the hedged portfolio is long the underlying and 
short the bank account. Conversely, for a short position in a vanilla put option, 
the hedged portfolio is short the underlying and long the bank account. Therefore, 
the call side of the butterfly spread rises when the borrowing rate is increased 
(cf. Figure[2]), and the put side rises when the borrowing rate is increased or stock 
borrowing fees are introduced (cf . Figures[2j[3]) . 



To see how well the penalty approximation (cf. Section[2| works in practice, we 
observe the difference between Problems |1.3| and |2.1| in dependence on the penalty 
parameter p, where we solve the former using policy iteration (see [101 H2] and 
Appendix[A| and the latter using Algorithm 3.2 We measure the difference in 



the maximum norm at time zero and keep the time and space grid fixed at M, 
N = 400. Clearly, Figure [4] shows that the observed convergence is of first order, 
which confirms the result of Corollary |2.6| (We find almost identical convergence 
rates when using (M, N) G {(600, 600), (1000, 1000), (900, 30), (30, 900)}, i.e. the 
observed penalty convergence is independent of the grid size.) 

Finally, we come to analyse the numerical efficiency of the penalty method devised 
in this paper. Especially in comparison with policy iteration, there are two relevant 
aspects to this: (i) the number of iterations needed to solve the nonlinear discrete 
system and (ii) the overall runtime of the method. 

The necessary numbers of iterations are listed in Table [2j and we see that, for 
the examples chosen, neither method ever needs more than four iterations; for the 
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Butterfly Spread — Unequal 
Borrowing/Lending Rates 




O 100 200 300 400 500 600 



Asset Price 



Figure 2. The value of a short position in a butterfly spread. When 
raising the borrowing rate, the value on the left hand side ( where a short 
position in the bank account is needed for hedging) increases and the right 
hand side does not change. When introducing a fee for stock borrowing, 
the right hand side (where a short position in the underlying is needed 
for hedging) increases and the left hand side does not change. 




Figure 3. The value of a short position in a butterfly spread. When 
decreasing the lending rate, the value on the right hand side (where a 
long position in the bank account is needed for hedging) increases and 
the left hand side does not change. When introducing a fee for stock 
borrowing, the right hand side (where a short position in the underlying 
is needed for hedging) increases again and the left hand side again does 
not change. 
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Policy Iteration 


n = 1 


n = 2 


n = 3 


n = 4 


M, N = 400 


90.50% 


9.50% 






M, TV = 1000 


91.20% 


8.80% 


- 


- 


M = 900, TV = 30 


99.67% 


0.33% 


- 


- 


M = 30, TV = 900 


- 


96.67% 


3.33% 


- 


Penalty Method (p = 4e03) 


n = 1 


n = 2 


n = 3 


n = 4 


M, TV = 400 






78.75% 


21.25% 


Af , AT = 1000 






78.40% 


21.60% 


M = 900, TV = 30 






92.56% 


7.44% 


M = 30, TV = 900 






66.67% 


33.33% 


Penalty Method (p = le06) 


n = 1 


n = 2 


n = 3 


n = 4 


M, TV = 400 






79.00% 


21.00% 


M, TV = 1000 






78.20% 


21.80% 


M = 900, iV = 30 






92.44% 


7.56% 


M = 30, TV = 900 






70.00% 


30.00% 



Table 2. 77ie table shows the number of iterations (denoted byn) needed 
when using policy iteration and Algorithm \3.2\ to solve the nonlinear dis- 
crete Problems \1.3\ and \2.1\ respectively. Since the number of required 
iterations might differ between time steps, the percentages denote how 
frequently a certain number occurred. Indisputably, both schemes con- 
verge astonishingly well (since generally n < 4), and, in many cases, 
policy iteration seems to even need only one step, whereas Algorithm \3.2\ 
predominantly needs three. 



Grid Size 


Policy 


Penalty (p = 4e03) 


Penalty (p = le06) 


M, TV = 400 


0.216s 


0.733s 


0.740s 


M , TV = 1000 


1.148s 


4.021s 


4.093s 


M = 900, TV = 30 


0.142s 


0.473s 


0.475s 


M = 30, TV = 900 


0.0624s 


0.116s 


0.121s 



Table 3. We see the computational time needed to solve Problems 1.3 
and \2.1\ using policy iteration and Algorithm \3.2\ respectively. The 
penalty scheme needs roughly three to four times longer than the policy 
iteration, a difference which reflects the different numbers of iterations 
listed in Table^ 



majority of cases, policy iteration needs only one iteration and Afgorithm |3.2| needs 
merely three. The overall runtimes of both approaches are listed in Table [3j and, in 
consistency with the number of iterations needed for solving the nonlinear system, 
we see that the penalty method runs up to three to four times as long as the policy 
iteration; nicely, the runtime of the penalty scheme does not depend on the size of 
the penalty parameter. 



Since we generally expect finite termination (as stated in Lemma 3.4 and Propo- 
sition 5.1|, one could, instead of proceeding as in Remark 5.2 also test for con- 
vergence by checking to what degree x n = x n+1 is satisfied. Interestingly, this is 
a disadvantage if a scheme converges in very few steps since one needs an extra 
step for verification once a certain accuracy is reached. (Using Remark 5.2 to test 



for convergence, Table [2] shows that, in most cases, the policy iteration reaches the 
desired accuracy in only one step. In [TD], where it is tested for \x" 
the authors generally need two steps.) 



< toL 
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Convergence in the Penalty Parameter 







Error 






Linear Regression 


0.02 






0.015 




'•^Slope = -0.9938 


0.01 






0.008 






0.006 






0.004 






0.003 







120 180 240 360 



600 
P 



Figure 4. We test the penalisation technique introduced in Section^ 
The plot (which is in log-log coordinates) shows the difference between the 
solutions to Problems ] 1.3\ and \2.1\ in dependence of the penalty parameter 
p > 0. As p is increased, the two problems converge nicely. Linear 
regression suggests the error is of order 0(l/p), confirming the analytical 



results ( cf. Theorem 2. 6 ) 



In total, for the current example, both approaches clearly work very well, with 
the policy iteration being computationally slightly more efficient (cf. Tablcs[2j[3j) 
than the penalty scheme. Since, analytically, both algorithms have very similar 



properties, i.e. both algorithms converge monotonically (cf. Lemma 3.3 [10 ) and 



in finitely many steps (cf. Lemma |3.4[ Proposition 5.1), it is hard to say generally 
which scheme should be faster. Since the frameworks of this paper and [TO] are not 
identicaQ (with the example of Section [5] fitting into both setups), the question of 
which algorithm is advantageous might simply be down to the exact nature of the 
problem. 

6. Conclusion 

We have shown that penalty methods can be a powerful tool for the solution of HJB 
equations with finitely many controls, and we have included a rigorous analysis of 
the convergence properties of the overall approach. More precisely, we have dis- 
cretised the HJB equation using standard finite differences and have approximated 
the resulting nonlinear discrete problem using a penalisation technique. If the dif- 
ferential operators of the equation obey a certain structure, we have been able to 
prove stability and convergence of the finite difference approximation in the viscos- 
ity sense and to estimate the penalisation error to be of order 0(1/ p), where p > 
is the penalty parameter. For the numerical solution of the penalised problem, we 
have devised an iterative method which has finite termination in theory and, for the 



^Loosely speaking, an HJB equation simply minimises over a (in our case finite) number of 
operators. The policy iteration presented in |10] needs a time derivative to appear in all operators. 
The penalty scheme of this paper can easily be extended to include time independent operators, 
like the identity function used when comparing the value of the solution to the payoff in the case 
of American-style options. 
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considered examples, converges in less than four steps in practice. There is a wide 
range of problems in mathematical finance matching our framework; in particular, 
we have presented numerical results for European option pricing in a model with 
different borrowing/lending rates and stock borrowing fees and have demonstrated 
the competitiveness of our approach. 



Appendix A. The Method of Policy Iteration 
We use this appendix to briefly describe the policy iteration scheme for the numer- 



ical solution of Problem 4.4 as devised in [TO]. Also, we show how convergence and 



finite termination can easily be deduced. 

In the situation of Problem |4.4[ we define 

1 



where I G K. denotes the identity matrix; in this new notation, Problem 4.4 
be reformulated: Find {V^)j e _\4 C l w such that 

yj-i _ yi 

(A.l) + min : (r, g)eS}= 0, j G M\{1}. 



Given some , j G the following algorithm solves equation (A.l I for 1 



Algorithm A.l. (Policy Iteration) Let x° G M. N be some starting value. Then, 
for known x n , n > 0, find x n+1 such that 

(A.2) (/+ k^f' a )x n+l = V\ 

where (r n ,q n ) G argmin{ 2l^| ,<T a; n : (r, q) G §}. 

Now, for n > 1, we have that 

(7 + fc</' CT )( a; "+ 1 -x") 



-1 n-1 



>(I + k% r B g q >°)x n+1 -(I + k 2l^ s ' q ' a )x n = 0, 

which, knowing that I + /cSl^ 9 ' a is an M-matrix, means that x n+1 > x n . Fur- 
thermore, since {(J + fcSlss' 7 ) -1 : ( r il) G S} is a finite set, we may deduce the 
convergence of (x n ) n& ^ and the existence of a A G N, independent of x°, such that 



x n+i _ x n j? or gjj n > A. If x n+1 = x n , expression (A.2) transforms into 

(I + kWgf' a )x n = V j 

& x n + k min { % r B q fx n : (r, q) G §} = V j , 
which means that x n solves ( AT] ). 



References 

[1] A. L. Amadori. Nonlinear integro-differential evolution problems arising in option pricing: A 
viscosity solution approach. Journal of Differential and Integral Equations, 16(7):787-811, 
2003. 

[2] G. Barles. Convergence of numerical schemes for degenerate parabolic equations arising in 

finance. In L. C. Rogers and D. Talay, editors, Numerical Methods in Finance, pages 1—21. 

Cambridge: Cambridge University Press, 1997. 
[3] G. Barles and E. R. Jakobsen. On the convergence rate of approximation schemes 

for Hamilton- Jacobi-Bellman equations. Mathematical Modelling and Numerical Analysis, 

36(l):33-54, 2002. 

[4] G. Barles and E. R. Jakobsen. Error bounds for monotone approximation schemes for para- 
bolic Hamilton- Jacobi-Bellman equations. Mathematics of Computation, 76(260):1861— 1893, 
2007. 



18 



J. H. WITTE AND C. REISINGER 



[5] A. Bcnsoussan and J. L. Lions. Applications of variational inequalities in stochastic control, 
volume 12 of Studies in Mathematics and its Applications. Amsterdam, New York, Oxford: 
North-Holland Pub. Co. , 1982. 

[6] Y. Z. Bergman. Option pricing with differential interest rates. The Review of Financial Stud- 
ies, 8(2):475-500, 1995. 

[7] S. Chaumont. Uniqueness to elliptic and parabolic Hamilton-Jacobi-Bcllman equations with 
non-smooth boundary. Coraptes Rendus Mathematique, 339(8) :555-560, 2004. 

[8] M. Fiedler. Special matrices and their applications in numerical mathematics. Lancaster: 
Nijhoff, 1986. 

[9] W. H. Fleming and H. M. Soner. Controlled Markov processes and viscosity solutions. New 

York: Springer, 2nd edition, 2005. 
[10] P. A. Forsyth and G. Labahn. Numerical methods for controlled Hamilton-Jacobi-Bcllman 

PDEs in finance. The Journal of Computational Finance, 11 (2): 1-44, 2007. 
[11] P. A. Forsyth and K. R. Vetzal. Quadratic convergence for valuing American options using a 

penalty method. SI AM Journal on Scientific Computing, 23(6):2095-2122, 2002. 
[12] P. A. Forsyth and K. R. Vetzal. Numerical methods for nonlinear PDEs in finance. 

Working paper, University of Waterloo, http://www.rmi.nus.edu.sg/csf/wcbpages/Authors 

/ scconddraft / forsyth-vetzal.pdf, 2010. 
[13] D. Duffie, N. Garleanu and L. H. Pedersen. Securities lending, shorting, and pricing. Journal 

of Financial Economics, 66(2-3) :307-339, 2002. 
[14] M. G. Crandall, H. Ishii and P.-L. Lions. User's guide to viscosity solutions of second order 

partial differential equations. Bulletin of the American Mathematical Society, 27(l):l-67, 

1992. 

[15] K. Ito and K. Kunisch. Semi-smooth Newton methods for variational inequalities of the first 
kind. Mathematical Modelling and Numerical Analysis, 37(l):41-62, 2003. 

[16] K. Ito and K. Kunisch. On a semi-smooth Newton method and its globalization. Mathematical 
Programming, 118(2):347-370, 2007. 

[17] H. J. Kushner. Numerical methods for stochastic control problems in continuous time. SIAM 
Journal on Control and Optimization, 28(5):99-1048, 1990. 

[18] H. J. Kushner and P. G. Dupuis. Numerical methods for stochastic control problems in con- 
tinuous time. New York: Springer, 2nd edition, 2001. 

[19] H. E. Lcland. Option pricing and replication with transactions costs. The Journal of Finance, 
40(5): 1283-13-1, 1985. 

[20] M. Avellaneda, A. Levy and A. Paras. Pricing and hedging derivative securities in markets 
with uncertain volatilities. Applied Mathematical Finance, 2(2):73-888, 1995. 

[21] M. S. Santos and J. Rust. Convergence properties of policy iteration. SIAM Journal on 
Control and Optimization, 42:2094-2115, 2004. 

[22] R. Seydel. Tools for computational finance. Univcrsitext. Berlin: Springer, 3rd edition, 2006. 

[23] Q. S. Song. Convergence of Markov chain approximation on generalized HJB equation and 
its applications. Automatica, 44(3):761-766, 2008. 

[24] R. S. Varga. Matrix iterative analysis. Berlin, London: Springer, 2000. 

[25] J. Wang and P. Forsyth. Maximal use of central differencing for Hamilton- Jacobi-Bellman 
PDEs in finance. SIAM Journal on Numerical Analysis, 46(3):1580-1601, 2009. 

[26] J. Yong and X. Y. Zhou. Stochastic controls: Hamiltonian systems and HJB equations. New 
York, London: Springer, 1999. 

Mathematical Institute, University of Oxford 
E-mail address: [witte , reisinge] (3maths.ox.ac.uk 



